The present invention relates to oil and gas production using computerized simulation of oil or gas reservoirs and of production facilities.
FIG. 1 schematically illustrates production facilities for extracting oil and/or gas from a reservoir 110. Reservoir 110 consists of subsurface rock formations below the earth surface 112. A wellbore 120.1 has been drilled into the reservoir and connected to surface facilities 130 a system of pipes and pumps 130 to pump the fluid (oil or gas, possibly mixed with water or mud) out of the reservoir. A casing has been inserted into wellbore 120.1, and perforations 134 were made in the casing to allow the wellbore to communicate with the reservoir. Another wellbore 120.2, connected to surface network 130 and having perforations 134, can optionally be used to inject water, steam, carbon dioxide or other chemicals into the reservoir to push the oil or gas towards wellbore 120.1 or make sticky oil more mobile.
Design of production facilities involves choice of wellbore locations and injected chemicals, well shapes and dimensions, pipe and pump networks, drilling schedules, production schedules, etc. The optimal choices of these and other parameters depend on physical properties of reservoir 110 such as pressures, saturations, porosities, and possibly others. Such choices are facilitated by computer simulation that simulates production for existing or hypothetical production facilities and production schedules. A suitable computer system 204 (FIG. 2) includes one or more computer processors 210 coupled to a memory 220. Typically, test measurements are made at selected locations in reservoir 110 or in existing production facilities to determine parameters of interest. The test data may include geological data from sources such as seismic analysis, rock cores, and well logs. These data may include data on rock porosities, permeabilities, and other information of interest, and allow determination of parameters of interest such as pressures or saturations at a number of locations. The test measurements are stored in memory 220 as shown at 230. A simulator program 240, stored in memory 220, is executed by processors 210 to read the test data 230 and provide simulated results which can be stored in memory 220 and/or displayed on a printer or computer monitor or other output device, or provided to another computer program or another computer over a network or otherwise, for use by oil field designers or other personnel.
A discrete model of reservoir 110 and the production facilities is created by computer 204 in memory 220 for use by simulation program 240 as illustrated in FIG. 3. More particularly, an imaginary grid of cells (grid blocks) 320 is superimposed on reservoir 110, with data indicating the locations of wellbores 120 and perforations 134. A node (not shown) is defined in each cell (typically, but not necessarily, at the center), and the pressure and possibly other values in the cell will be modeled as the pressure or other values at the node. The wellbore portions above the perforations are shown as 120.1W, 120.2W and will be called just “wells” herein. Nodes 324 are defined along each well 120W (i.e. 120.1W or 120.2W) in the computer model; these nodes represent point locations at which the pressures or other values will modeled. A well may have one or more nodes 324. Connections 330 represent the well portions between the nodes 324. Similar nodes and connections (not shown) can represent the surface network 130. Other nodes 325 may be defined to represent point locations in the perforated sections of the wellbores (below the wells). Connections 331 represent the wellbore portions between the nodes 325 and between the top wellbore node 325 and the first well node 324 (the “bottom-hole” node). A flow of material between a perforated cell 320 (i.e. a cell perforated by a wellbore) and the wellbore will be modeled as the flow between the perforated cell and a node 325, or a flow between the cell and the corresponding bottom-hole node 324.
A coordinate system (XYZ) is chosen to represent various locations. A time axis 304 with time points t1, t2, t3 . . . is defined in computer 204.
A set of equations is chosen to model the oil or gas field, including hydraulics equations between nodes in the wells, wellbores, and the surface network; Darcy and non-Darcy flow between cells 320; perforation equations between cells 320 and nodes 325; fluid component mass balance equations for nodes and cells; constraint equations, and possibly others. See U.S. Pat. No. 7,668,707 issued Feb. 23, 2010 to Watts, III et al., incorporated herein by reference. The equations may include differential equations, and the differential equations are discretized using known techniques. The equations can be symbolically written as:F(Vari,j,k,l)=0  (1)
where Vari,j,k,l represents pressures, flow rates, or other variables at the nodes and connections at different times. The subscript “i” represents the time point ti on axis 304, and the subscripts j, k, l represent X, Y, and Z values respectively. The values Vari,j,k,l are treated as independent variables.
Since some physical parameters can be expressed through others, there is some latitude as to which parameters are chosen as the independent variables (primary variables) in equations (1). Further, the equations (1) depend on the discretization method, and in particular on whether any particular Var value is measured at a cell 320 or a node 324 or 325, or a segment 330 or 331. Such choices will affect how efficiently the computer can solve the equations. The term “efficiently” relates to minimal use of the computer resources, and particularly the processor and memory. It is desirable to reduce the execution time and memory amount needed to solve the equations while providing a solution that serves as a good approximation to the real world.
There can be millions of different locations and time points, so the equation (1) may represent millions or billions of scalar equations. In a typical procedure, the equations (1) are solved for each of consecutive times ti, i.e. for time t1, then for time t2, and so on. The solution for time t1 depends on the test measurements stored as shown at 230 in FIG. 2. The solution for each subsequent time ti is found using the previously determined Var values for time ti−1 and possibly earlier times and possibly using the test measurements 230.
Equations (1) are typically non-linear, and for each time ti these equations can be solved approximately using the Newton-Raphson method illustrated in FIG. 4. This is an iterative method in which the non-linear equations (1) are replaced by linear equations at each iteration. The successive iterations use different linear equations to provide hopefully better and better approximate solutions at each iteration. To simplify the notation, let us re-write the equations (1) as:F(x)=0
Suppose xn is a current iterative solution in the Newton-Raphson method. Then the next iteration xn+1 is obtained as follows. The Jacobian matrix JF(xn) is calculated; this is the matrix of partial derivatives of F(x) evaluated at the current iteration xn. Then the following linear equation is solved for xn+1:JF(xn)(xn+1−xn)=−F(xn)
This equation can be alternatively written as:JF(xn)(xn+1)=JF(xn)(xn)−F(xn)  (2)orF′L(xn+1)=b  (3)
where F′L=JF, and b is the constant equal to the right-hand side of equation (3).
This process is illustrated in FIG. 4. The symbol Var.NT denotes the current approximation xn. Since ti is fixed, the Vari,j,k,l will sometimes be written as simply Varj,k,l. The values Var.NT of the last iteration are taken as the approximate solution Vari,j,k,l of the equations (1) for the time ti.
The initial iterated values Var.NT is chosen at step 410, and this can be any desired values, e.g. the values Vari−1,j,k,l found in the previous iteration i−1 if i>1, or to some values obtained from the initial measurements 230 if i=1.
At step 414, the matrix JF(Var.NT) and the valueb=JF(Var.NT)Var.NT−F(Var.NT)
are determined for the equations (2)-(3). At step 418, these equations are solved, i.e. the following equation is solved:F′L(Varj,k,l)=b  (4)
Let us denote the solution as Var*j,k,l. If this is the last iteration, then the solution Var*j,k,l is taken as the solution of (1), i.e. the method of FIG. 4 terminates for ti. See step 422. If other iterations are desired, then the solution Var*j,k,l is used as the next value of Var.NT for the next iteration (step 426), and steps 414-418 are repeated. Any number of iterations can be used to obtain better and better approximations of the solution of (1). The iterations can be terminated when F(Var*j,k,l) is close to zero in some sense.
Equations (4) may include billions of equations with billions of individual variables shown collectively as Varj,k,l. Therefore, computer system 240 solves equations (4) approximately. An exemplary approximate method is a Schwarz method illustrated in FIG. 5. For simplicity, the function F′L will be denoted simply as f, and Varj,k,l as x. Thus, equation (4) can be written as:f(x)=b  (4′)
The function f represents many scalar linear functions, and the variable x represents just as many scalar components. The Schwarz method will be illustrated on the example in which each of x and f has 1000 scalar components, denoted as x1, x2, . . . , x1000 for x, and denoted as f1, f2, . . . , f1000 for f.
At step 610, the values x1, . . . , x1000 are divided into two overlapping groups of values (more than two groups can also be used). In the example discussed below, the two groups are:G1: x1, . . . , x600G2: x400, . . . x1000
Also, the scalar functions f1, . . . , f1000 are divided into two overlapping groups, with the same number of functions in each group as there are values in G1, G2. For example, the two groups of functions can be:G1′: f1, . . . , f600G2′: f400, . . . , f1000
We will use the notation G1′(x) to denote the values f1(x) through f600(x); the notation G2′(x) will denote the values f400(x) through f1000(x).
The Schwarz method is an iterative method, and the iterative solution will be denoted as x.sc. The value x.sc output by the Schwarz method is the solution x.
At step 620, the initial iterative value x.sc is determined using any suitable technique. For example, x.sc can be set to the value of Var.NT at step 418 of the current Newton-Raphson iteration of FIG. 4.
At steps 630-640, the next iterated value is determined for x.sc as follows.
At step 630, the x values from G1\G2, i.e. the values x601 through x1000, are set to the corresponding values of the current iterative value x.sc, and the following equation is solved for x1 through x600:G1′(x)=b1  (5)
where b1 is the “b” portion corresponding to G1′. In other words, the following equation is solved:G1′(x1, . . . , x600, x.sc601, . . . , x.sc1000)=b1  (5′)
The solution will be denoted x1*, . . . , x600*.
The values x1* through x399* will be used as the corresponding values of the next iterative value of x.sc as explained below.
At step 640, the following equation is solved for x400* through x1000*:G2′(x)=b2  (6)
where b2 is the b portion corresponding to G2. In this equation, the values x1, . . . , x399 are constant. In the additive Schwarz method, these values are set to the corresponding values of the current iterative value x.sc. In the multiplicative Schwarz method, these values are set to the corresponding values x1*, . . . , x399* found at step 630. The solution of (6) will be denoted x600*, . . . , x1000*. The values x1* through x1000* will be collectively denoted as x*.
If this is the last iteration, then x*can be taken as the solution of (4′), i.e. as the value Var*j,k,l at step 418. See step 650. If other iterations are desired (step 660), then the solution x* is used as the next iterated value of x.sc, and steps 630-640 are repeated. Any number of iterations can be used to obtain better and better approximations of the solution of (4′).
In some variations of the Schwarz method, iterations other than the first iteration are solved for the residual b−f(x) rather than x. The residual is updated at each iteration. In particular, the overlap values x400*, . . . , x600* found at step 630 are used to update the residual for step 640.
In a variation, the components x400*, . . . , x600* of x* are obtained for step 650 by combining the corresponding values obtained at steps 630 and 640, i.e. by combining the corresponding components of solutions of equations (5) and (6) so as to minimize the residual's norm ∥b−f(x*)∥. A known method for combining these components is GMRES (Generalized Minimum Residual method).
Clearly, the computational complexity of solving the equations (1) depends on the choice of variables Var and the groups G1, G2, G1′, G2′. Generally, the reservoirs are quite different from wells with respect to the physics models, so the reservoirs and the wells are usually in different groups. More particularly, each group G1, G2 can be defined through cells 320 or nodes 324 or 325: each group consists of variables pertaining to a group of cells 320 or nodes 324 or 325 or the related segments 330 or 331, and each function group G1′, G2′ consists of the equations for the corresponding cells and nodes. Thus, the reservoir cells can be placed in group G1 (i.e. in the node group defining G1), the wells in group G2, and the perforation cells 320 (the cells with perforations) and wellbore nodes 325 can be the overlap between G1 and G2. This scheme is used in NEXUS® Reservoir Simulation Software system available from Halliburton of Houston, Tex. In NEXUS®, the variables Var are pressures and component compositions at nodes 324 and 325, and pressures and component masses in the cells 320.
Background information can also be found in the following publications incorporated herein by reference (no representation is being made as to whether or not these publications are prior art relative to the present disclosure):                U.S. patent application publication 2009/0294122 A1 (Dec. 3, 2009), inventors Hansen et al.        U.S. Pat. No. 7,668,707 B2 (Feb. 23, 2010), inventors Watts, III et al.        U.S. Pat. No. 7,516,056 B2 (Apr. 7, 2009), inventors Wallis et al.        U.S. patent application publication 2010/0286971 A1 (Nov. 11, 2010), inventors Middya et al.        UK patent application GB 2 434 235 A (Jul. 18, 2007), applicant Schlumberger Technology Corporation.        European patent application EP 2 192 507 A1 (Jun. 2, 2010), applicant Maersk Olie OG Gas A/S.        PCT application WO 2006/020952 A2 (Feb. 23, 2006), applicants Saudi Arabian Oil Company and ARAMCO Services Company.        H. Kazemi, “Numerical Simulation of Water-Oil Flow in Naturally Fractured Reservoirs”, Society of Petroleum Engineers Journal, December 1976, pages 317-326.        Rafael Bru et al., “Overlapping additive and multiplicative Schwarz iterations for H-matrices”, Linear Algebra and its Applications, Elsevier, 2003.        